Robust Detection of Dynamic Community Structure in Networks 
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We describe techniques for the robust detection of community structure in some classes of time- 
dependent networks. Specifically, we consider the use of statistical null models for facilitating the 
principled identification of structural modules in semi-decomposable systems. Null models play 
an important role both in the optimization of quality functions such as modularity and in the 
subsequent assessment of the statistical validity of identified community structure. We examine the 
sensitivity of such methods to model parameters and show how comparisons to null models can help 
identify system scales. By considering a large number of optimizations, we quantify the variance of 
network diagnostics over optimizations ('optimization variance') and over randomizations of network 
structure ('randomization variance'). Because the modularity quality function typically has a large 
number of nearly- degenerate local optima for networks constructed using real data, we develop a 
method to construct representative partitions that uses a null model to correct for statistical noise 
in sets of partitions. To illustrate our results, we employ ensembles of time-dependent networks 
extracted from both nonlinear oscillators and empirical neuroscience data. 



PACS numbers: 89.75.Fb, 89.75.Hc, 87.19.L-, 87.18.Vf 

Many social, physical, technological, and bio- 
logical systems can be modeled as networks com- 
posed of numerous interacting parts [lj. As an 
increasing amount of time-resolved data has be- 
come available, it has become increasingly impor- 
tant to develop methods to quantify and charac- 
terize dynamic properties of temporal networks 
[2]. Generalizing the study of static networks, 
which are typically represented using graphs, 
to temporal networks entails the consideration 
of nodes (representing entities) and/or edges 
(representing ties between entities) that vary in 
time. As one considers data with more compli- 
cated structures, the appropriate network anal- 
yses must become increasingly nuanced. In the 
present paper, we discuss methods for algorithmic 
detection of dense clusters of nodes (i.e., commu- 
nities) by optimizing quality functions on mul- 
tilayer network representations of temporal net- 
works [3, 4j. We emphasize the development and 
analysis of different types of null-model networks, 
whose appropriateness depends on the structure 
of the networks one is studying as well as the con- 
struction of representative partitions that take 
advantage of a multilayer network framework. To 



illustrate our ideas, we use ensembles of time- 
dependent networks from the human brain and 
human behavior. 



INTRODUCTION 

Myriad systems have components whose interactions 
(or the components themselves) change as a function of 
time. Many of these systems can be investigated using 
the framework of temporal networks, which consist of sets 
of nodes and/or edges that vary in time [2|. The formal- 
ism of temporal networks is convenient for studying data 
drawn from areas such as person-to-person communica- 
tion (e.g., via mobile phones [5j|6]), one-to-many infor- 
mation dissemination (such as Twitter networks [7 ) , cell 
biology, distributed computing, infrastructure networks, 
neural and brain networks, and ecological networks [2]. 
Important phenomena that can be studied in this frame- 
work include network constraints on gang and criminal 
activity [HI [9] , political processes [10l [11] , human brain 
function [4j [12] , human behavior [13] , and financial struc- 
tures [HJCE5]. 

Time-dependent complex systems can have densely 
connected components in the form of cohesive groups of 



2 




1 1 1 » 

Time 

FIG. 1. (Color Online) An important property of many real- 
world networks is community structure, in which there exist 
cohesive groups of nodes such that a network has stronger con- 
nections within such groups than it does between such groups. 
Community structure often changes in time, which can lead 
to the rearrangement of cohesive groups, the formation of new 
groups, and the breakup of existing groups. 

nodes known as 'communities' (see Fig. [I]), which can be 
related to a system's functional modules [16 , 17 . A wide 
variety of clustering techniques have been developed to 
identify communities, and they have yielded insights in 
the study of the committee structure in the United States 
Congress [18 j, functional groups in protein interaction 
networks [15] , functional modules in brain networks [4], 
and more. A particularly successful technique for iden- 
tifying communities in networks p~6j [20| is optimization 
of a quality function known as 'modularity' [21], which 
recently has been generalized for detecting communities 
in time-dependent and multiplex networks [3] . 

Modularity optimization allows one to algorithmically 
partition a network's nodes into communities such that 
the total connection strength within groups of the parti- 
tion is more than would be expected in some null model. 
However, modularity optimization always yields a net- 
work partition (into a set of communities) as an output 
whether or not a given network truly contains modu- 
lar structure. Therefore, application of subsequent diag- 
nostics to a network partition is potentially meaningless 
without some comparison to benchmark or null-model 
networks. That is, it is important to establish whether 
the partition(s) obtained appear to represent meaning- 
ful community structures within the network data or 
whether they might have reasonably arisen at random. 
Moreover, robust assessment of network organization de- 
pends fundamentally on the development of statistical 
techniques to compare structures in a network derived 
from real data to those in appropriate models (see, e.g., 
Ref. [22 ). Indeed, as the constraints in null models and 
network benchmarks become more stringent, it can be- 
come possible to make stronger claims when interpreting 
organizational structures such as community structure. 

In the present paper, we examine null models in time- 
dependent networks and investigate their use in the al- 
gorithmic detection of cohesive, dynamic communities in 
such networks (see Fig. Indeed, community detection 
in temporal networks necessitates the development of null 
models that are appropriate for such networks. Such null 
models can help provide bases of comparison at various 



stages of the community-detection process, and they can 
thereby facilitate the principled identification of dynamic 
structure in networks. Indeed, the importance of devel- 
oping null models extends beyond community detection, 
as such models make it possible to obtain statistically 
significant estimates of network diagnostics. 

Our dynamic network null models fall into two cat- 
egories: optimization null models, which we use in 
the identification of community structure; and post- 
optimization null models, which we use to examine the 
identified community structure. We describe how these 
null models can be selected in a manner appropriate to 
known features of a network's construction, identify po- 
tentially interesting network scales by determining values 
of interest for structural and temporal resolution param- 
eters, and inform the choice of representative partitions 
of a network into communities. 



METHODS 
Community Detection 

Community-detection algorithms provide ways to de- 
compose a network into dense groups of nodes called 
'modules' or 'communities'. Intuitively, a community 
consists of a set of nodes that are connected among 
one another more densely than they are to nodes in 
other communities. A popular way to identify commu- 
nity structure is to optimize a quality function, which 
can be used to measure the relative densities of intra- 
community connections versus inter-community connec- 
tions. See [16] [20] [23] for recent reviews on network com- 
munity structure and [24H27] for discussions of various 
caveats that should be considered when optimizing qual- 
ity functions to detect communities. 

One begins with a network of N nodes and a given 
set of connections between those nodes. In the usual 
case of single-layer networks (e.g., static networks with 
only one type of edge), one represents a network using 
an TV x TV adjacency matrix A. The element of the 
adjacency matrix indicates a direct connection or 'edge' 
from node i to node j, and its value indicates the weight 
of that connection. The quality of a hard partition of 
A into communities (whereby each node is assigned to 
exactly one community) can be quantified using a quality 
function. The most popular choice is modularity [16] [20] 

|2T][2E1[25] 

= Yl\ Ai i - 7 p ij] s (9i,9j) , (!) 

ij 

where node i is assigned to community g^ node j is as- 
signed to community gj, the Kronecker delta S(gi,gj) = 1 
if gi = gj and it equals otherwise, 7 is a resolution 
parameter (which we will call a structural resolution pa- 
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FIG. 2. (Color Online) Methodological considerations im- 
portant in the investigation of dynamic community structure 
in temporal networks. (A) Depending on the system under 
study, a single network layer (which is represented using an 
ordinary adjacency matrix with an extra index to indicate the 
layer) might by definition only allow edges from some subset 
of the complete set of node pairs, as is the case in the depicted 
chain-like graph. We call such a situation partial connectivity. 
(B) Although the most common optimization null model em- 
ploys random graphs (e.g., the Newman-Girvan null model, 
which is closely related to the configuration model [TJ [16]), 
other models can also provide important insights into net- 
work community structure. (C) After determining a set of 
partitions that maximize the modularity Q (or a similar qual- 
ity function), it is interesting to test whether the community 
structure is different from, for example, what would be ex- 
pected with a scrambling of time layers (i.e., a temporal null 
model) or node identities (i.e., a nodal null model) [3]. 



rameter) , and P^ is the expected weight of the edge con- 
necting node i to node j under a specified null model. 
The choice 7 = 1 is very common, but it is important to 
consider multiple values of 7 to examine groups at multi- 
ple scales [I6j [30j [31] . Maximization of Qo yields a hard 
partition of a network into communities such that the 



total edge weight inside of modules is as large as possible 
(relative to the null model and subject to the limitations 
of the employed computational heuristics, as optimizing 
Qo is NP-hard [II [2UJ [32] ) . 

Recently, the null model in the quality function (??) 
has been generalized so that one can consider sets of L 
adjacency matrices, which are combined to form a rank-3 
adjacency tensor A that can be used to represent time- 
dependent or multiplex networks. One can thereby define 
a multilayer modularity (also called 'multislice modular- 
ity') [3] ' 
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{(Aiji - liPiji) Sir + S ij uj j i r }S(g i i 1 g jr ) 



(2) 

where the adjacency matrix of layer I has components 
Aiji, the element P^\ gives the components of the corre- 
sponding layer-/ matrix for the optimization null model, 
7/ is the structural resolution parameter of layer Z, the 
quantity gu gives the community assignment of node i in 
layer Z, the quantity gj r gives the community assignment 
of node j in layer r, the element U0ji r gives the connection 
strength (i.e., an 'interlayer coupling parameter', which 
one can call a temporal resolution parameter if one is us- 
ing the adjacency tensor to represent a time-dependent 
network) from node j in layer r to node j in layer Z, the 
total edge weight in the network is /1 = \ J2j r fyr, the 
strength (i.e., weighted degree) of node j in layer I is 
Kji = kji + Cj/, the intra-layer strength of node j in layer 
I is kji — "^ Ji Am, and the inter-layer strength of node j 
in layer I is c jt = J2 r u jir- 

Equivalent representations that use other notation can, 
of course, be useful. For example, multilayer modularity 
can be recast as a set of rank-2 matrices describing con- 
nections between the set of all nodes across layers [e.g., 
for spectral partitioning [29, 33j[34]]. One can similarly 
generalize Q for higher-rank tensors, which one would use 
when studying community structure in networks that are 
both time-dependent and multiplex, through appropriate 
specification of inter-layer coupling tensors. 



Network Diagnostics 

To characterize multilayer community structure, we 
compute four example diagnostics for each hard parti- 
tion: the modularity Q, the number of modules n, the 
mean community size s (which is equal to the number of 
nodes in the community and is proportional to 1/n), and 
the stationarity ( [35 . To compute we calculate the 
autocorrelation function U(t,t + m) of two states of the 
same community G(t) at m time steps (i.e., m network 
layers) apart: 



U(t,t + m) 



= \G(t)nG(t + m)\ 
= |G(*)UG(* + m)| 



(3) 
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where \G(t) n G(t + m)\ is the number of nodes that are 
members of both G(t) and G(£+ra), and |G(t)UG(£+ra)| 
is the number of nodes in the union of the community at 
times t and t + m. Defining to to be the first time step 
in which the community exists and t' to be the last time 
in which it exists, the stationarity of a community is [35] 

. Et>(M + l) (4) 
s t'-t " w 

This gives the mean autocorrelation over consecutive 
time steps [36] . 

In addition to these diagnostics, which are defined us- 
ing the entire multilayer community structure, we also 
compute two example diagnostics on the community 
structures of the component layers: the mean single- 
layer modularity (Q s ) and the variance var(Q s ) of the 
single- layer modularity over all layers. The single- layer 
modularity Q s is defined as the static modularity qual- 
ity function, Q s = Y.ij[ A ij ~ lPij]S(9i,9j), computed 
for the partition g that we obtained via optimization of 
the multilayer modularity function Q. We have chosen 
to use a few simple ways to help characterize the time 
series for Q S1 though of course other diagnostics should 
also be informative. 



DATA SETS 

We illustrate the concept and uses of dynamic net- 
work null models using two example network ensembles: 

(1) 75-time-layer brain networks drawn from each of 20 
human subjects and (2) behavioral networks with about 
150 time layers drawn from each of 22 human subjects. 
Importantly, the use of network ensembles makes it pos- 
sible to examine robust structure (and also its variance) 
over multiple network instantiations. We have previously 
examined both data sets in the context of neuroscientific 
questions [4] [13]. In this paper, we use them as illus- 
trative examples for the consideration of methodological 
issues in the detection of dynamic communities in tem- 
poral networks. 

These two data sets, which provide examples of differ- 
ent types of network data, illustrate a variety of issues 
in network construction: (1) node and edge definitions, 

(2) complete versus partial connectivity, (3) ordered ver- 
sus categorical nodes, and (4) confidence in edge weights. 
In many fields, determining the definition of nodes and 
edges is itself an active area of investigation [37 . See, for 
example, several recent papers that address such ques- 
tions in the context of large-scale human brain networks 
[38H43] and in networks more generally [44 . Another im- 
portant issue is whether to examine a given adjacency 
matrix in an exploratory manner or to impose struc- 
ture on it based on a priori knowledge. For example, 
when nodes are categorical, one might represent their re- 



lations using a fully connected network and then iden- 
tify communities of any group of nodes. However, when 
nodes are ordered — and particularly when they are in 
a chain of weighted nearest-neighbor connections — one 
expects communities to group neighboring nodes in se- 
quence, as typical community-detection methods are un- 
likely to yield many out-of-sequence jumps in community 
assignment. The issue of confidence in the estimation of 
edge weights is also very important, as it can prompt an 
investigator to delete edges from a network when their 
statistical validity is questionable. A closely related is- 
sue is how to deal with known or expected missing data, 
which can affect either the presence or absence of nodes 
themselves or the weights of edges [45H48] . 

Data Set 1: Brain Networks 

Our first data set contains categorical nodes with par- 
tial connectivity and variable confidence in edge weights. 
The nodes remain unchanged in time, and edge weights 
are based on covariance of node properties. This covari- 
ance structure is non-local in the sense that weights exist 
between both topologically neighboring nodes and topo- 
logically distant nodes [49] [50] . This property has been 
linked in other dynamical systems to behaviors such as 
chimera states, which coherent and incoherent regions 
coexist [5TH53] , Another interesting feature of this data 
set is that it is drawn from an experimental measurement 
with high spatial resolution (on the order of centimeters) 
but relatively poor temporal resolution (on the order of 
seconds). 

As described in more detail in Ref. [4] , we construct an 
ensemble of networks (20 individuals over 3 experiments, 
which yields 60 multilayer networks) that represent the 
functional connectivity between large regions of the hu- 
man brain. In these networks, N = 112 centimeter-scale, 
anatomically distinct brain regions are our (categorical) 
network nodes. We study the temporal interaction of 
these nodes — such interactions are thought to underly 
cognitive function — by first measuring their activity ev- 
ery 2 seconds during simple finger movements using func- 
tional magnetic resonance imaging (fMRI). We cut these 
regional time series into time slices (which yield layers 
in the multilayer network) of roughly 3-minute duration. 
Each such layer corresponds to a time series whose length 
is 80 units. 

To estimate the interactions (i.e., edge weights) be- 
tween nodes, we calculate a measure of statistical sim- 
ilarity between regional activity profiles [54]. Using a 
wavelet transform, we extract frequency-specific activity 
from each time series in the range 0.06-0.12 Hz. For each 
time layer I and each pair of regions i and j, we define the 
weight of an edge connecting region i to region j using 
the coherence between the wavelet-coefficient time series 
in each region, and these weights form the elements of 
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a weighted, undirected temporal network W with com- 
ponents Wiji = Wju. The magnitude- squared coherence 
Gij between time series i and j is a function of frequency. 
It is defined by the equation 



Gij(f) 



\F i3 (f)f 



(5) 



where Fu(f) and Fjj(f) are the power spectral density 
functions of i and j, respectively, and Fij(f) is the cross- 
power spectral density function of i and j. We let Hij 
denote the mean of Gij(f) over the frequency band of 
interest, and the weight of edge Wiji is equal to com- 
puted for layer I. 

We use a false-discovery rate [55] to threshold connec- 
tions whose coherence values are not significantly greater 
than that expected at random. This yields a multilayer 
network A with components Aiji (i.e., a rank-3 adjacency 
tensor). The nonzero entries in retain their weights. 
We couple the layers of A^i to one another with tempo- 
ral resolution parameters of weight ujji r between node j 
in layer r and node j in layer I. In this paper, we let 
Ujir = uj G [0.1, 40] be identical between each node j in a 
given layer with itself in nearest-neighbor layers. (In all 
other cases, ujji r = 0.) 

In Fig. we show an example time layer from A^ji 
for a single subject in this experimental data. In this 
example, the statistical threshold is evinced by the set of 
matrix elements set to 0. Because brain network nodes 
are categorical, one can apply community detection algo- 
rithms in these situations to identify communities com- 
posed of any set of nodes. (Note that the same node 
from different times can be assigned to the same com- 
munity even if the node is assigned to other communities 
at intervening times.) One biological interpretation of 
network communities in brain networks is that they rep- 
resent groups of nodes that serve distinct cognitive func- 
tions (e.g., vision, memory, etc.) that can vary in time 

pa Eg. 



Data Set 2: Behavioral Networks 

Our second data set contains ordered nodes that re- 
main unchanged in time. The network topology in this 
case is highly constrained, as edges are only present be- 
tween consecutive nodes. (We call this 'nearest-neighbor' 
coupling.) Another interesting feature of this data set is 
that the number of time slices is an order of magnitude 
larger than the number of nodes in a slice. 

As described in more detail in Ref. [13 , we construct 
an ensemble of 66 behavioral networks from 22 individu- 
als and 3 experimental conditions. These networks repre- 
sent a set of finger movements in the same simple motor 
learning experiment from which we constructed the brain 
networks in data set 1. Subjects were instructed to press 



a sequence of buttons corresponding to a sequence of 12 
pseudo-musical notes shown to them on a screen. 

Each node represents an interval between consecutive 
button presses. A single network layer consists of N = 11 
nodes (i.e., there is one interval between each pair of 
notes), which are connected in a chain via weighted, undi- 
rected edges. In Ref. [T3] , we examined the phenomenon 
of motor 'chunking', which is a fascinating but poorly- 
understood phenomenon in which groups of movements 
are made with similar inter- movement durations. (This 
is similar to remembering a phone number in groups of 
a few digits or grouping notes together as one masters 
how to play a song.) For each experimental trial I and 
each pair of inter- movement intervals i and j, we de- 
fine the weight of an edge connecting inter- movement i 
to inter- movement j as the normalized similarity in inter- 
movement durations. The normalized similarity between 
nodes i and j is defined as 



Piji 



d\ — dij 



ijl 



(6) 



where d^\ is the absolute value of the difference of lengths 
of the zth and jth inter- movement time intervals in trial 
/ and di is the maximum value of diji in trial I. These 
weights yield the elements Wiji of a weighted, undirected 
multilayer network W. Because finger movements occur 
in series, inter-movement i is connected in time to inter- 
movement i ± 1 but not to any other inter-movements 
i + n for \n\ ^ 1. 

To encode this conceptual relationship as a network, 
we set all non-contiguous connections in Wiji to and 
thereby construct a weighted, undirected chain network 
A^i. In Fig.[3j3, we show an example trial layer from A^i 
for a single subject in this experimental data. We couple 
layers of A^i to one another with weight u)ji r , which gives 
the connection strength between node j in experimental 
trial r and node j in trial /. In a given instantiation of the 
network, we again let U0ji r = lj G [0.1,40] be identical for 
all nodes j for all connections between nearest-neighbor 
layers. (Again, Uji r = in all other cases.) Because inter- 
movement nodes are ordered, one can apply community- 
detection algorithms to identify communities of nodes in 
sequence. Each community represents a motor 'chunk'. 



RESULTS 
Modularity-Optimization Null Models 

After constructing a multilayer network A with el- 
ements Aiji, it is necessary to select an optimization 
null model P^i in equation (??). The most common 
modularity-optimization null model used in undirected, 
single-layer networks is the Newman-Girvan null model 



6 



Brain Network Layer g Behavioral Network Layer 
lol 




Coherence 



Dynamic Communities 
C in Brain Networks 



Similarity 
Dynamic Communities 
D in Behavioral Networks 




Communities 



Communities 



FIG. 3. (Color Online) Network layers and community as- 
signments from two example data sets: (A) a brain network 
based on correlations between blood-oxygen-level-dependent 
(BOLD) signals [4] and (B) a behavioral network based on 
similarities in movement times during a simple motor learn- 
ing experiment [13 . We use these data sets to illustrate situ- 
ations with categorical nodes and ordered nodes, respectively. 
In the bottom panels, we show community assignments ob- 
tained using multilayer community detection for (C) the brain 
networks and (D) the behavioral networks. 
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where ki = ■ is the strength of node i and m = 
\ Aij. The definition (??) can be extended to mul- 



tilayer networks using 



ijl 



kukji 
2mi 



(8) 



where ku = ^ Mji is the strength of node i in layer / and 
mi = | Y^ij Aiji. Optimization of Q using the null model 
(??) identifies partitions of a network into groups that 
have more connections (in the case of binary networks) 
or higher connection densities (in the case of weighted 
networks) than would be expected for the distribution of 
connections (or connection densities) expected in a null 
model. We use the notation A 1 for the layer-/ adjacency 
matrix composed of elements A^\ and the notation to 
denote the layer-/ null- model matrix with elements Pij\. 
See Fig. [SjA. for an example layer A 1 from a multilayer 
behavioral network and Fig. [4^3 for an example instanti- 
ation of the Newman-Girvan null model P l . 



Optimization Null Models for Ordered Node Networks 

The Newman-Girvan null model is particularly useful 
for networks with categorical nodes, in which a connec- 
tion between any pair of nodes can occur in theory. How- 
ever, when using a chain network of ordered nodes, it is 
useful to consider alternative null models. For example, 
in an ordinary network (i.e., one that is represented using 
an adjacency matrix), one can define 



P 



(9) 



where p is the mean edge weight of the chain network 
and A[- is the binarized version of A^ , in which nonzero 
elements of Aij are set to 1 and zero-valued elements 
remain unaltered. Such a null model can also be defined 
for multilayer networks: 



Piji = PiA't 



ijl •> 



(10) 



where pi is the mean edge weight in layer / and A! { - x is 
the binarized version of A^i. The optimization of Q using 
this null model identifies partitions of a network whose 
communities have a larger strength than the mean. See 
Fig. for an example of this chain null model for 
the behavioral network layer shown in Fig. |4]A.. 

In Fig. [4)3, we illustrate the effect that the choice of 
optimization null model has on the modularity values Q 
of the behavioral networks as a function of the structural 
resolution parameter. (Throughout the manuscript, we 
use a Louvain-like locally-greedy algorithm to maximize 
the multilayer modularity quality function [571158].) The 
Newman-Girvan null model gives decreasing values of Q 
for 7 G [0.1,2.1], whereas the chain null model produces 
lower values of Q, which behaves in a qualitatively differ- 
ent manner for 7 < 1 versus 7 > 1. To help understand 
this feature, we plot the number and mean size of com- 
munities as a function of 7 in Figs. |4^1 and [4^. As 7 is 
increased, the Newman-Girvan null model yields network 
partitions that contain progressively more communities 
(with progressively smaller mean size). The number of 
communities that we obtain in partitions using the chain 
null model also increases with 7, but it does so less grad- 
ually. For 7 <C 1, one obtains a network partition con- 
sisting of a single community of size Ni = 11; for 7 ^> 1, 
each node is instead placed in its own community. For 
7 = 1, nodes are assigned to several communities whose 
constituents vary with time (see, for example, Fig. J3p). 

The above results highlight the sensitivity of network 
diagnostics such as Q, n, and s to the choice of an op- 
timization null model. It is important to consider this 
type of sensitivity in the light of other known issues, such 
as the extreme near-degeneracy of quality functions like 
modularity [23] . Importantly, the use of the chain null 
models provides a clear delineation of network behavior 
in this example into three regimes as a function of 7: 
a single community with variable Q (low 7), a variable 
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FIG. 4. (Color Online) Modularity-optimization null models. (A) Example layer A 1 from a behavioral network. (B) Newman- 
Girvan and (C) chain null models for the layer shown in panel (A). (D) Optimized multilayer modularity value Q, (E) 
number of communities n, and (F) mean community size s for the complete multilayer behavioral network employing the 
Newman-Girvan (black) and chain (red) optimization null models as a function of the structural resolution parameter 7. (G) 
Optimized modularity value Q, (H) number of communities n, and (/) mean community size s for the multilayer behavioral 
network employing chain optimization null models as a function of the effective fraction £ m i(7) of edges that have larger weights 
than their null-model counterparts. We averaged the values of Q, n, and s over the 3 different 12-note sequences and C = 100 
optimizations. Box plots in (D-F) indicate quartiles and 95% confidence intervals over the 22 individuals in the study. The 
error bars in panels (G-I) indicate a standard deviation from the mean. In some instances, this is smaller than the line width. 
The temporal resolution-parameter value is uj — 1. 



number of communities as Q reaches a minimum value 
(7 « 1), and a set of singleton communities with min- 
imum Q (high 7). This illustrates that it is crucial to 
consider a null model appropriate for a given network, 
as it can provide more interpret able results than just us- 
ing the usual choices (such as the Newman-Girvan null 
model). 

The structural resolution parameter 7 can be trans- 
formed so that it measures the effective fraction of edges 
£(7) that have larger weights than their null- model coun- 
terparts [31] . One can define a generalization of £ to mul- 
tilayer networks, which allows one to examine the behav- 
ior of the chain null model near 7 = 1 in more detail. 
For each layer /, we define a matrix X z (7) with elements 
Xiji(j) = Aiji — jPiji, and we then define c x (7) to be 
the number of elements of X z (7) that are less than 0. 
We sum c x (7) over layers in the multilayer network to 
construct ^(7). The transformed structural resolution 



parameter is then given by 



6nl(7) = ~x 



ml 



(7) 



(Amin) 



(A, 



(Amin) 



(11) 



where A m i n is the value of 7 for which the network still 
forms a single community in the multilayer optimization 
and A max is the value of 7 for which the network still 
forms N singleton communities in the multilayer opti- 
mization. (We use Roman typeface in the subscripts in 
c* Y and £ m i to emphasize that we are describing multi- 
layer objects and, in particular, that the subscripts do not 
represent indices.) In Figs.[4]3-I, we report the optimized 
(i.e., maximized) modularity value, the number of com- 
munities, and the mean community size as functions of 
the transformed structural resolution parameter £ m i(7)- 
(Compare these plots to Figs. [Ip-F.) For all three diag- 
nostics, the apparent transition points seem to be more 
gradual as a function of £ m i(7) than they are as a function 
of 7. For systems like the present one that do not exhibit 
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a pronounced, nontrivial plateau in these diagnostics as 
a function of a structural resolution parameter, it might 
be helpful to have a priori knowledge about the expected 
number or sizes of communities (see, e.g., [13 ) to help 
guide further investigation. 



Optimization Null Models for Networks Derived from Time 
Series 



Although the Newman- Girvan null model can be used 
in networks with categorical nodes, such as the brain 
networks in data set 1 (see Fig.[5]A.), it does not take ad- 
vantage of the fact that these networks are derived from 
similarities in time series. Accordingly, we generate sur- 
rogate data to construct two dynamic network null mod- 
els for community detection that might be particularly 
appropriate for networks derived from time-series data. 

First, we note that a simple null model (which we call 
'Random') for time series is to randomize the elements 
of the time-series vector for each node before computing 
the similarity matrix (see Fig. [5}3) [59]. However, the 
resulting time series do not have the mean or variance of 
the original time series, and this yields a correlation- or 
coherence-based network with very low edge weights. To 
preserve the mean, variance, and autocorrelation func- 
tion of the original time series, we employ a surrogate- 
data generation method that scrambles the phase of time 
series in Fourier space [60]. Specifically, we assume that 
the linear properties of the time series are specified by 
the squared amplitudes of the discrete Fourier transform 



\S(u)\' 



1 



v-i 

E i27ruv/V 
s v e 

v o 



(12) 



where s v denotes an element in a time series of length V. 
(That is, V is the number of elements in the time-series 
vector.) We construct surrogate data by multiplying the 
Fourier transform by phases chosen uniformly at random 
and transforming back to the time domain: 



1 y_1 



2-Kkv/V 



(13) 



where a u G [0, 2tt) are chosen independently and uni- 
formly at random. [61] This method, which we call the 
Fourier transform (FT) surrogate (see Fig.[5]C), has been 
used previously to construct covariance matrices [62] and 
to characterize networks [63]. A modification of this 
method, which we call the amplitude- adjusted Fourier 
transform (A AFT) surrogate, allows one to also retain 
the amplitude distribution of the original signal [64] (see 
Fig. [5J3). One can alter nonlinear relationships between 
time series while preserving linear relationships between 
time series by applying an identical shuffling to both time 



series; one can alter both linear and nonlinear relation- 
ships between time series by applying independent shuf- 
flings to each time series [60 . 

We demonstrate in Fig. [5^1 that, among the four null 
models that we consider, the mean coherence of pairs of 
FT surrogate series match most closely to that of the 
original data. Pairs of Random time series have the 
smallest mean coherence, and pairs of A AFT surrogate 
series have the next smallest. The fact that the A AFT 
surrogate is less like the real data (in terms of mean co- 
herence) than the simpler FT surrogate might stem from 
a rescaling step [62] that causes the power spectrum to 
whiten (i.e., the step flattens the power spectral density) 
[65] . In Figs.^-H, we show three diagnostics (optimized 
modularity, mean community size, and number of com- 
munities) as a function of the structural resolution pa- 
rameter 7 for the various optimization null models. We 
note that the Newman- Girvan null model produces the 
smallest Q value and a middling community size, whereas 
the surrogate time series models produce higher Q val- 
ues and more communities of smaller mean size. The 
Random null model produces the largest value of Q and 
the fewest communities, which is consistent with the fact 
that it contains the smallest amount of shared informa- 
tion (i.e., mean coherence) with the real network. 



Post-Optimization Null Models 

After identifying the partition (s) that maximize mod- 
ularity, one might wish to determine whether the iden- 
tified community structure is significantly different from 
that expected under other null hypotheses. For example, 
one might wish to know whether any temporal evolu- 
tion is evident in the dynamic community structure (see 
Fig. |2p). To do this, one can employ post- optimization 
null models, in which a multilayer network is scrambled in 
some way to produce a new multilayer network. One can 
then maximize the modularity of the new network and 
compare the resulting community structure to that ob- 
tained using the original network. Unsurprisingly, one's 
choice of post-optimization null model should be influ- 
enced by the question of interest, and it should also be 
constrained by properties of the network under examina- 
tion. We explore such influences and constraints using 
our example networks. 



Intra-Layer and Inter-Layer Null Models 

There are various ways to construct connectional null 
models (i.e., intra-layer null models), which randomize 
the connectivity structure within a network layer (A 1 ) 
[4] [66]. For binary networks, one can obtain ensembles 
of random graphs with the same mean degree as that of 
a real network using Erdos-Renyi random graphs [1 , and 
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FIG. 5. (Color Online) Modularity-optimization null models for time series. (^4) Example coherence matrix A 1 averaged over 
layers from a brain network. (B) Random time shuffle, (C) Fourier transform (FT) surrogate, and (D) amplitude-adjusted 
Fourier transform (AAFT) surrogate null models P l averaged over layers. (E) Coherence of each matrix type averaged over 
subjects, scans, and layers. We note that the apparent lack of structure in (B) is partially related to its significantly decreased 
coherence in comparison to the other models. (F) Optimized modularity values Q, (G) number of communities n, and (H) 
mean community size s for the multilayer brain network employing the Newman-Girvan (black), random time-shuffle (blue), 
FT surrogate (gray), and AAFT surrogate (red) optimization null models as functions of the structural resolution parameter 7. 
We averaged the values of these diagnostics over 3 different scanning sessions and C = 100 optimizations. Box plots indicate 
quartiles and 95% confidence intervals over the 20 individuals in the study. The temporal resolution parameter is uj — 1. 



ensembles of weighted random networks can similarly be 
constructed from weighted random graph models [67] . To 
retain both the mean and distribution of edge weights, 
one can employ a permutation-based connectional null 
model that randomly rewires network edges with no ad- 
ditional constraints by reassigning uniformly at random 
the entire set of matrix elements Aiji in the Ith layer 
(i.e., the matrix A*). Other viable connectional null mod- 
els include ones that preserve degree [2TJ [68] or strength 
[69] distributions, or — for networks based on time-series 
data — preserve length, frequency content, and ampli- 
tude distribution of the original time series [70]. In this 
section, we present results for a few null models that are 
applicable to a variety of temporal networks. We note, 
however, that this is a fruitful area of further investiga- 
tion. 

We employ two connectional null models specific for 
the broad classes of networks represented by the brain 
and behavioral networks that we use as examples in this 
paper. The brain networks provide an example of time- 
dependent similarity networks, which are weighted and 
either fully connected or almost fully connected [3T] . 
(The brain networks have some entries in their cor- 
responding adjacency tensors because we have removed 
edges with weights that are not statistically significant 
[4] .) We therefore employ a constrained null model that is 
constructed by randomly rewiring edges while maintain- 
ing the empirical degree distribution [68]. In Fig.[6]Al, we 



demonstrate the use of this null model to assess dynamic 
community structure. Importantly, this constrained null 
model can be used in principle for any binary or weighted 
network, though it does not take advantage of specific 
structure (aside from strength distribution) that one 
might want to exploit. For example, the behavioral net- 
works have chain-like topologies, and it is desirable to 
develop models that are specifically appropriate for such 
situations. (One can obviously make the same argument 
for other specific topologies.) We therefore introduce a 
highly constrained connectional null model that is con- 
structed by reassigning edge weights uniformly at random 
to existing edges. This does not change the underlying 
binary topology. (That is, we preserve network topology 
but scramble network geometry.) We demonstrate the 
use of this null model in Fig. [6j31. 

In addition to intra-layer null models, one can also em- 
ploy inter-layer null models — such as ones that scramble 
time or node identities [4 . For example, we construct a 
temporal null model by randomly permuting the order of 
the network layers. This temporal null model can be used 
to probe the existence of significant temporal evolution 
of community structure. One can also construct a nodal 
null model by randomly permuting the inter-layer edges 
that connect nodes in one layer to nodes in another. Af- 
ter the permutation is applied, an inter-layer edge can, 
for example, connect node i in layer t with node j 7^ i in 
layer t + 1 rather than being constrained to connect each 
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node i in layer t with itself in layer t + 1. One can use 
this null model to probe the importance of node iden- 
tity in network organization. We demonstrate the use 
of our temporal null model in row 2 of Fig. [6j and we 
demonstrate the use of our nodal null model in row 3 of 
Fig.© 

Calculation of Diagnostics on Real Versus Null-Model 
Networks 

We characterize the effects of post-optimization null 
models using four diagnostics: maximized modularity Q, 
the number of communities n, the mean community size 
s, and the stationarity ( (see the section titled 'Network 
Diagnostics' for definitions). Due to the possibly large 
number of partitions with nearly optimal Q [24] , the val- 
ues of such diagnostics vary over realizations of a compu- 
tational heuristic for both the real and null-model net- 
works. (We call this optimization variance.) The null- 
model networks also have a second source of variance 
(which we call randomization variance) from the myriad 
possible network configurations that can be constructed 
from a randomization procedure. We note that a third 
type of variance — ensemble variance — can also be 
present in systems containing multiple networks. In the 
example data sets that we discuss, this represents vari- 
ability among experimental subjects. 

We test for statistical differences between the real 
and null- model networks as follows. We first compute 
C = 100 optimizations of the modularity quality function 
for a network constructed from real data and then com- 
pute the mean of each of the four diagnostics over these 
C samples. This yields representative values of the di- 
agnostics. We then maximize modularity for C different 
randomizations of a given null model (i.e., 1 optimization 
per randomization) and then compute the mean of each 
of the four diagnostics over these C samples. For both of 
our example data sets, we perform this two-step proce- 
dure for each network in the ensemble (60 brain networks 
and 66 behavioral networks; see 'Methods'). We then in- 
vestigate whether the set of representative diagnostics 
for the networks constructed from real data are different 
from those of appropriate ensembles of null-model net- 
works. To address this issue, we subtract the diagnostic 
value for the null model from that of the real network 
for each subject and experimental session. We then use 
one-sample t-tests to determine whether the resulting dis- 
tribution differs significantly from 0. We show our results 
in Fig. [6j 

Results depend on all three factors (the data set, the 
null model, and the diagnostic), but there do seem to 
be some general patterns. For example, the real net- 
works exhibit the most consistent differences from the 
nodal null model for all diagnostics and both data sets 
(see row 2 of Fig. |6|. For both data sets, the variance 



of single-layer modularity in the real networks is con- 
sistently greater than those for all three null models, 
irrespective of the mean (see the final two columns of 
Figs. [6]A.,B); this is a potential indication of the statis- 
tical significance of the temporal evolution. However, 
although optimized modularity is higher in the real net- 
work for both data sets, the number of communities is 
higher in the set of brain networks and lower in the set 
of behavioral networks. Similarly, in comparison to the 
connectional null model, higher modularity is associated 
with a smaller mean community size in the brain net- 
works but a larger mean size in the behavioral networks 
(see row 1 of Fig. [6|. These results demonstrate that 
the three post-optimization null models provide different 
information about the network structure of the two sys- 
tems and thereby underscores the critical need for further 
investigations of null-model construction. 

Structural and Temporal Resolution Parameters 

When optimizing multilayer modularity, we must 
choose (or otherwise derive) values for the structural res- 
olution parameter 7 and the temporal resolution param- 
eters uj. By varying 7, one can tune the size of commu- 
nities within a given layer: large values of 7 yield more 
communities, and small values yield fewer communities. 
A systematic method for how to determine values of ujji r 
has not yet been discussed in the literature. In princi- 
ple, one could choose different ujji r values for different 
nodes, but we focus on the simplest scenario in which 
the value of ujji r = uj is identical for all nodes j and all 
contiguous pairs of layers I and r (and is otherwise 0). 
In this framework, the temporal resolution parameter uj 
provides a means of tuning the number of communities 
discovered across layers: high values of uj yield fewer com- 
munities, and low values yield more communities. It is 
beneficial to study a range of parameter values to exam- 
ine the breadth of structural (i.e., intra-layer [24| [25 | f7Tj ) 
and temporal (i.e., inter-layer) resolutions of community 
structure, and some papers have begun to make progress 
in this direction [3j II [[3j HQ [72] . 

To characterize community structure as a function of 
resolution-parameter values (and hence of system scales), 
we quantify the quality of partitions using the mean value 
of optimized Q. To do this, we examine the constitution 
of the partitions using the mean similarity over C op- 
timizations, and we compute partition similarities using 
the z-score of the Rand coefficient [73] . For comparing 
two partitions a and /3, we calculate the Rand z-score in 
terms of the network's total number of pairs of nodes M, 
the number of pairs M a that are in the same community 
in partition a, the number of pairs Mp that are in the 
same community in partition /?, and the number of pairs 
w a f3 that are assigned to the same community both in 
partition a and in partition j3. The z-score of the Rand 
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FIG. 6. (Color Online) Post-optimization null models. We compare four multilayer diagnostics (optimized modularity, number 
of communities, mean community size, and stationarity) and two single-layer diagnostics (mean and variance of Q s ) for (A) 
brain and (B) behavioral networks with the connectional (row 1), nodal (row 2), and temporal (row 3) null-model networks. 
Box plots indicate quartiles and 95% confidence intervals over the individuals and experimental conditions. The structural 
resolution parameter is 7 = 1 and the temporal resolution parameter is uj — 1. 



coefficient comparing these two partitions is 



1 



M 



(14) 



where <J Waf3 is the standard deviation of w a p (as in [73]). 
Let the mean partition similarity z denote the mean value 
of z a p over all possible partition pairs for a ^ f3. 

In Fig.[7| we show both z and optimized Q as a function 
of 7 and uj in both brain and behavioral networks. The 
highest modularity values occur for low 7 and high uj. 
The mean partition similarity is high for large 7 in the 
brain networks, and it is high for both small and large 7 in 
the behavioral networks. Interestingly, in both systems, 
the partition similarity when 7 = uj = 1 is lower than it is 
elsewhere in the (7, uj) parameter plane, so the variability 
in partitions tends to be large at this point. Indeed, as 
shown in the second row of Fig. [7| modularity exhibits 
significant variability for 7 = uj = 1 compared to other 
resolution-parameter values. 

It is useful to be able to determine the ranges of 7 
and uj that produce community structure that is signif- 
icantly different from a particular null model. One can 
thereby use null models to probe resolution-parameter 
values at which a network displays interesting structures. 
This could be especially useful in systems for which one 
wishes to identify length scales (such as a characteris- 
tic mean community size) or time scales [4j |35] [74j [75] 
directly from data. 



In Fig. [8j we show examples of how the difference 
between diagnostic values for real and null-model net- 
works varies as a function of 7 and uj. As illustrated in 
panels (A) and (B), the brain and behavioral networks 
both exhibit a distinctly higher mean optimized modu- 
larity than the associated nodal null-model network for 
7 « uj « 1. Interestingly, this roundly peaked differ- 
ence in Q is not evident in comparisons of the real net- 
works to temporal null-model networks (see Figs. [8p,D), 
so resolution-parameter values (and hence system scales) 
of potential interest might be more identifiable by com- 
parison to nodal than to temporal null models in these 
examples. It is possible, however, that defining tempo- 
ral layers over a longer or shorter duration would yield 
identifiable peaks in the difference in Q. 

The differences in the Rand z-score landscapes are 
more difficult to interpret, as the values of mean parti- 
tion similarity z are much larger in the real networks for 
some resolution-parameter values (positive differences; 
red) but are much larger in the null-model networks for 
other resolution-parameter values (negative differences; 
blue). The clearest situation occurs when comparing 
the brain's real and temporal null-model networks (see 
Fig. [8p), as the network built from real data exhibits 
a much larger value of z (and hence much more con- 
sistent optimization solutions) than the temporal null- 
model networks for high values of 7 (i.e., when there 
many communities) and low uj (i.e., when there is weak 
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FIG. 7. (Color Online) Optimized modularity Q and Rand z-score as functions of the resolution parameters 7 and uj for the 
(A) brain and (B) behavioral networks. The top row shows the mean value of maximized Q over C = 100 optimizations and 
the mean partition similarity z over all possible pairs of the C partitions. The bottom row shows the variance of maximized 
Q over the optimizations and the variance of the partition similarity over all possible pairs of partitions. The results shown in 
this figure come from a single individual and experimental scan, but we obtain qualitatively similar results for other individuals 
and scans. Note that the axis scalings are nonlinear. 



temporal coupling) . These results are consistent with the 
fact that weak temporal coupling in a multilayer network 
facilitates greater temporal variability in network parti- 
tions across time. Such variability appears to be signif- 
icantly different than the noise induced by scrambling 
time layers. These results suggest potential resolution 
values of interest for the brain system, as partitions are 
very consistent across many optimizations. For example, 
it would be interesting to investigate community struc- 
ture in these networks for high 7 (e.g., 7 ~ 40) and low 
uj (e.g., uj rH 0.1). At these resolution values, one can 
identify smaller communities with greater temporal vari- 
ability than the communities identified for the case of 
7 = uj = 1 [4]. 

The optimization and randomization variances appear 
to be similar in the brain and behavioral networks (see 
rows 2-3 in every panel of Fig. |8| not only in terms of 
their mean values but also in their distribution in the 
part of the (7, uj) parameter plane that we examined. In 
particular, the variance in Q is larger in the real net- 
works precisely where the mean is also larger, so mean 
and variance are likely either dependent on one another 
or on some common source. Importantly, such depen- 
dence influences the ability to draw statistical conclusions 
because it is possible that the points in the (7, uj) plane 
with the largest differences in mean are not necessarily 
the points with the most significant differences in mean. 

We also find that the dependencies of the diagnostics 
on 7 and uj are consistent across subjects and scans, sug- 
gesting that our results are ensemble-specific rather than 
individual-specific. 



Examination of Data Generated from a Dynamical 
System 

Real- world data is often clouded by unknown or math- 
ematically undefinable sources of variance, so it is also 
important to examine data sets generated from dynami- 
cal systems (or other models). Because we are concerned 
with time-dependent networks, we consider an example 
consisting of time-dependent data generated by a well- 
known dynamical system. 

We construct a network of Kuramoto oscillators, in 
which the phase 6i (t) of the i th oscillator evolves in time 
according to 

-jjT=Ui + J2 KAijSiniOj -Ot), i e {1, . . . , N} , (15) 

3 

where uji is the natural frequency of oscillator z, the ma- 
trix A gives the binary coupling between each pair of os- 
cillators, and k is a positive real constant that indicates 
the strength of the coupling. We draw the frequencies uji 
from a Gaussian distribution with mean and standard 
deviation 1. In our simulations, we use a time step of 
r = 0.1, a constant of k — 0.2, and a network size of 
N = 128. 

Kuramoto oscillators have been studied in the context 
of various network topologies and geometries [51 53] [76]- 
[78] and from both the component and ensemble perspec- 
tives [79 . We are interested in networks with dynamic 
community structure. Following Refs. [77] [80], we im- 
pose a well-defined community structure in which each 
community is composed of 16 nodes. In each time step, 
each node has 13 connections with nodes in its own com- 
munity and 1 connection with nodes outside of its com- 
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FIG. 8. (Color Online) Differences, as a function of 7 and cj, between the real networks and the (A,B) nodal and (C,D) temporal 
null models for maximized modularity Q and partition similarity z for the (A,C) brain and (B,D) behavioral networks. The 
first row in each panel gives the difference in the mean values of the diagnostic variables between the real and null-model 
networks. Panels (A,B) show the results for Q — Q n and z — z n , and panels (C,D) show the results for Q — Q 1 and z — z t . The 
quantities Q and z again denote the modularity and partition similarity of the real network, Q n and z n denote the modularity 
and partition similarity of the nodal null-model network, and Q* and z t denote the modularity and partition similarity of the 
temporal null-model network. The second row in each panel gives the difference between the optimization variance of the real 
network and the randomization variance of the null-model network for the same diagnostic variable pairs. The third row in 
each panel gives the difference in the optimization variance of the real network and the optimization variance of the null-model 
network for the same diagnostic variable pairs. We show results for a single individual and scan in the experiment, but results 
are qualitatively similar for other individuals and scans. Note that the axis scalings are nonlinear. 



munity (see Fig. |9]A.). 

To quantify the temporal evolution of synchronization 
patterns, we define a set of temporal networks given by 
the time-dependent correlation between pairs of oscilla- 



tors: 

<p ij (t) = (\cos[9 i (t)-9 j (t)]\) , (16) 

where the angular brackets indicate an average over 20 
simulations. As time evolves from time step t = to 
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FIG. 9. (Color Online) Dynamic community detection in a 
network of Kuramoto oscillators. (A, top) The coupling ma- 
trix between N = 128 phase oscillators contains 8 communi- 
ties, each of which has 16 nodes. (A, bottom) Over time, oscil- 
lators synchronize with one another. Color indicates the mean 
phase correlation between oscillators, where hotter (darker 
gray) colors indicate stronger correlations. (B) Phase corre- 
lation between oscillators as a function of time. The mean 
phase correlation between oscillators in the same community 
(dashed red curve) increases faster than the mean phase corre- 
lation between all oscillators in the system (solid gray curve). 
Regime I encompasses the first 50 time steps, and regime II 
emcompasses the subsequent 50 time steps. (C) Variance of 
maximized multilayer modularity (left), number of communi- 
ties (middle), and partition similarity z (right) over 100 op- 
timizations of the multilayer modularity quality function for 
the temporal network in regime II as a function of the struc- 
tural resolution parameter for uj — 1. The shaded gray area 
indicates values of the structural resolution parameter that 
provide variance in the number of communities. (D) Exam- 
ple partition of the temporal network in regime II at 7 = 1.5, 
which occurs near the troughs in panel (C). (E) Example 
partition of the temporal network in regime I at 7 = 1.5. 
(F) Number of communities as a function of time for (left) 
the temporal network in regime I and (right) its correspond- 
ing temporal null model. We averaged values over C = 100 
optimizations of multilayer modularity. 



t = 100, oscillators tend to synchronize with other oscil- 
lators in their same community more quickly than with 
oscillators in other communities (see Fig. [9^3). 

To examine the performance of our multilayer 
community-detection techniques in this example, we 
compute Aiji = Aijt = 4>i,j(t) and using the multi- 



layer extension of the Newman- Girvan null model P^i 
given in Eq.( ??). We separately optimize Q for two 
temporal regimes: (1) regime I (with t G {1, ...,50}), 
for which synchronization within communities increases 
rapidly; and (2) regime II (with t G {51, . . . , 100}), for 
which within-community synchronization level is roughly 
constant but global synchronization still increases gradu- 
ally. We set uo = 1 and probe the effects of the structural 
resolution parameter 7 in regime II. In Figs. [9p,D, we 
illustrate that one can identify the value of 7 that best 
uncovers the underlying hard-wired connectivity using 
troughs in the optimization variance of several diagnos- 
tics (e.g., maximized modularity, number of communi- 
ties, and mean partition similarity). 

We probe the community structure in regime I us- 
ing the value of 7 that best uncovered the underlying 
hard-wired connectivity in regime II. We observe tempo- 
ral changes of community structure at early time points, 
as evidenced by the large number of communities for 
t G {1,...,5} (see Figs. |9^1,F). Importantly, the tem- 
poral dependence of community number on t is not ex- 
pected from a post-optimization temporal null model (see 
the right panel of Fig. [9]F). We obtain qualitatively sim- 
ilar results when we optimize the multilayer modularity 
quality function over the entire temporal network with- 
out separating the data into two regimes. 

Our results illustrate that one can use dynamic com- 
munity detection to uncover the resolution of inherent 
hard-wired structure in a data set extracted from the 
temporal evolution of a dynamical system and that post- 
optimization null models can be used to identify regimes 
of unexpected temporal dependence in network structure. 



Dealing With Degeneracy: Constructing 
Representative Partitions 

The multilayer modularity quality function has numer- 
ous near-degeneracies, so it is important to perform many 
instantiations when using a non-deterministic computa- 
tional heuristic to optimize modularity [24]. In doing 
this, an important issue is how (and whether) to distill 
a single representative partition from a (possibly very 
large) set of C partitions [81 j. In Fig. 



10, we illustrate 



a new method for constructing a representative partition 
based on statistical testing in comparison to null models. 

Consider C partitions from a single layer of an example 
multilayer brain network (see Fig. 10 A.). We construct a 



nodal association matrix T, where the element indi- 
cates how many times nodes i and j have been assigned to 
the same community (see Fig.fTofi)). We then construct a 
null-model association matrix T r based on random per- 



mutations of the original partitions (see Fig. 10 3). That 
is, for each of the C partitions, we reassign nodes uni- 
formly at random to the n communities of mean size s 
that are present in the selected partition. For every pair 
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of nodes i and j, we let TT- be the number of times these 
two nodes have been assigned to the same community 
in this permuted situation (see Fig. 



10 



5). The values T[- 
then form a distribution for the expected number of times 
two nodes are assigned to the same partition. Using an 
example with C = 100, we observe that two nodes can be 
assigned to the same community up to about 30 times out 
of the C partitions purely by chance. To be conservative, 
we remove such 'noise' from the original nodal associa- 
tion matrix T by setting any element whose value is 
less than the maximum entry of the random association 
matrix to (see Fig. [Top). This yields the thresholded 
matrix T', which retains statistically significant relation- 
ships between nodes. 

We use a Louvain-like algorithm to perform C op- 
timizations of the single-layer modularity Qo for the 
thresholded matrix T' . Interestingly, this procedure typ- 
ically extracts identical partitions for each of these opti- 
mizations in our examples (see Fig. [Top). This method 
therefore allows one to deal with the inherent near- 
degeneracy of the modularity quality function and pro- 
vides a robust, representative partition of the original 
example brain network layer (see Fig. [Top). 

We apply the same method to multilayer networks 



(see Fig. 11) to find a representative partition of (1) a 
real network over C optimizations, (2) a temporal null- 
model network over C randomizations, and (3) a nodal 
null-model network over C randomizations. Using these 
examples, we have successfully uncovered representative 
partitions when they appear to exist (e.g., in the real net- 
works and the temporal null-model networks) and have 
not been able to uncover a representative partition when 
one does not appear to exist (e.g., in the nodal null- 
model network, for which each of the 112 brain nodes 
is placed in its own community in the representative par- 
tition). We also note that the representative partitions in 
the temporal null-model and real networks largely match 
the original data in terms of both sizes and number of 
communities. These results indicate the potential of this 
method to uncover meaningful representative partitions 
over optimizations or randomizations in multilayer net- 
works. 



CONCLUSIONS 

In this paper, we discussed methodological issues in the 
determination and interpretation of dynamic community 
structure in multilayer networks. We also analyzed the 
behavior of several null models used for optimizing qual- 
ity functions (such as modularity) in multilayer networks. 

We described the construction of networks and the ef- 
fects that certain choices can have on the use of both 
optimization and post-optimization null models. We in- 
troduce novel modularity-optimization null models for 
two cases: (1) networks composed of ordered nodes 
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FIG. 10. (Color Online) Constructing representative parti- 
tions for an example brain network layer. (A) Partitions ex- 
tracted during C optimizations of the quality function Q. (B) 
The N x N nodal association matrix T, whose elements in- 
dicate the number of times node i and node j have been as- 
signed to the same community. (C) The N xN random nodal 
association matrix T r , whose elements indicate the number 
of times node i and node j are expected to be assigned to 
the same community by chance. (D) The thresholded nodal 
association matrix T', where elements with values less than 
those expected by chance have been set to 0. (E) Partitions 
extracted during C = 100 optimizations of the single-layer 
modularity quality function Q s for the matrix T from panel 
(D). Note that each of the C optimizations yields the same 
partition. (F) Visualization of the representative partition 
given in (E) [82]. We have reordered the nodes in the ma- 
trices in panels (A-E) for better visualization of community 
structure. 



(a 'chain null model') and (2) networks constructed 
from time-series similarities (FT and AAFT surrogates) . 
We studied 'connectional', 'temporal', and 'nodal' post- 
optimization null models using several multilayer diag- 
nostics (optimized modularity, number of communities, 
mean community size, and stationarity) as well as novel 
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FIG. 11. (Color Online) Representative partitions of multi- 
layer brain networks for an example subject and scan. (A) 
Partitions extracted for C — 100 optimizations of the qual- 
ity function Q on the real multilayer network (112 nodes x 
25 time windows, which yields 2800 nodes in total). Parti- 
tions extracted for C randomizations for the (B) temporal 
and (C) nodal null-model networks. (D) Partitions extracted 
for C optimizations of the quality function Q of the thresh- 
olded nodal association matrix for the (D) real, (E) temporal 
null-model, and (F) nodel null-model networks. Note that 
the partitioning is robust to multiple optimizations. We have 
reordered the nodes in each column for better visualization of 
community structure. The structural resolution parameter is 
7 = 1, and the temporal resolution parameter is uj = 1. 



single-layer diagnostics (in the form of measures based on 
time series for optimized modularity). To investigate the 
utility of such considerations for model-generated data, 
we also applied our methodology to time-series data gen- 
erated from coupled Kuramoto oscillators. 

We examined the dependence of optimized modular- 
ity and partition similarity on structural and temporal 
resolution parameters as well as the influence of their 
variances on putative statistical conclusions. Finally, we 
described a simple method to address the issue of near- 
degeneracy in the landscapes of modularity and similar 
quality functions using a method to construct a robust, 
representative partition from a network layer. 

The present paper illustrates what we feel are impor- 
tant considerations in the detection of dynamic commu- 
nities. As one considers data with increasingly compli- 



cated structures, network analyses of such data must be- 
come increasingly nuanced, and the purpose of this paper 
has been to discuss and provide some potential starting 
points for how to try to address some of these nuances. 
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